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Abstract 



The Method of Fundamental Solutions (MFS) is a popular tool to solve Laplace and Helmholtz boundary value 
problems. Its main drawback is that it often leads to ill-conditioned systems of equations. In this paper we investigate 
for the interior Helmholtz problem on analytic domains how the singularities (charge points) of the MFS basis 
functions have to be chosen such that approximate solutions can be represented by the MFS basis in a numerically 
stable way. For Helmholtz problems on the unit disc we give a full analysis which includes the high frequency (short 
wavelength) limit. For more difficult and nonconvex domains such as crescents we demonstrate how the right choice 
of charge points is connected to how far into the complex plane the solution of the boundary value problem can be 
analytically continued, which in turn depends on both domain shape and boundary data. Using this we develop a 
recipe for locating charge points which allows us to reach error norms of typically 10 -11 on a wide variety of analytic 
domains. At high frequencies of order only 3 points per wavelength are needed, which compares very favorably to 
boundary integral methods. 
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1. Introduction 

The Method of Fundamental Solutions (MFS), also known as the charge simulation method or the method 
of auxiliary sources, is a well known method for solving Laplace or Helmholtz boundary value problems 
(BVPs). The idea is to approximate the solution by fundamental solutions of the Laplace or Helmholtz 
equation whose singularities lie outside the domain. Consider the boundary value problem 



where fl C K 2 = C is a simply connected planar domain with analytic boundary dfl. Recall that the solution 
is unique if and only if k 2 is not a Dirichlet eigenvalue (of the Laplace operator) for the domain; physically 
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waves 



Au + k 2 u = in fi, 
u = v on dfl, 
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this is a resonance effect. The idea of the MFS is to approximate u by a linear combination of fundamental 
solutions of the form 

. N 

U (x)««W(x) = l2a^(A|x-y J .|), Yj e R 2 \U, (2) 

where -ffg 1 ^ is a Hankel function of the first kind of order zero, and N is the number of approximating functions 
each of which is associated with a charge point y •. It is well known that satisfies the Helmholtz equation 
in C\{0} with a singularity at zero. It is common to choose charge points lying on a smooth curve; we then 
find it enlightening to interpret the MFS as a discretization of the single layer potential representation of u 
as follows. Let T be a closed curve enclosing ft such that dist(T, 90) := min{|x — y|, x e dfl, y e T} > 0, 
then given a density g G i x (r) we may write 

u(x) w - [ H^\k\x-s\)g(s) ds, xe!l. (3) 

If g(s) = X^jLi a j<K s — Yj) f° r somc point set {y^} € I\ where 5 is the Dirac delta, we recover the MFS 
formulation (2). Note that the irregular Bessel function Y , or Hankel may be used instead of in 

the MFS [7,10]; see Remark 1 below. 

An overview about the history of this method and its applications is given in [8] . The rate of convergence 
of the MFS for the Laplace BVP was investigated in [3,12,13,14,15]. It turns out that if the boundary data is 
analytic one can achieve exponential convergence for the MFS for the Laplace problem on analytic domains 
if the charge points y ■ are suitably chosen. 

One of the main drawbacks of the method is that in the end systems of equations or linear least squares 
problems have to be solved that are often ill-conditioned. The effects of this ill-conditioning on the quality 
of the solution have been investigated for the Laplace problem in [17,18]. In this paper we investigate 
more closely under what conditions on the points y 3 a numerically stable representation of an approximate 
solution of the Helmholtz problem (1) as linear combination of fundamental solutions is possible. It turns 
out that this depends on how far into the complex plane a solution of (1) can be analytically continued. The 
importance of this in the context of scattering problems has already been observed [20,21] (and references 
in [21]). 

Our work also has consequences for the numerical solution of more challenging and widely-applicable 
PDE problems that are closely related to the one we study. We have in mind i) finding eigenmodes of the 
Laplace operator in £1 with homogeneous boundary conditions (where the MFS has been used at low [7] and 
very high eigenvalue [2]), and ii) scattering of time-harmonic waves (the exterior Helmholtz boundary value 
problem in R 2 \ Q). In both these situations the boundary data is almost always analytic: in problem i) it 
is homogeneous and in ii) a plane wave or point source. 

We will study convergence of the MFS approximation in the boundary error norm 

t=\\u {N) -v\\ L 2 (m) . (4) 
By applying [19, Eq. 7], this controls the interior error of the solution as follows, 

\\u^-u\\ L2m <^|| U W- W || i2(8n) , (5) 

where d := min.,- \k 2 — Ej\/Ej, the domain's Dirichlet eigenvalues are Ej, and Cn is a domain-dependent 
constant. This shows that for any fixed nonresonant k, we may use the boundary norm. 

In Section 2 we give rigorous results for the convergence and the numerical stability of the MFS for 
Helmholtz problems on the unit disc, with analytic boundary data, using charge points on a concentric 
circle. We then present a heuristic model for behavior in finite-precision arithmetic and show it explains well 
numerical results observed at both low and high wavenumbers. A key conclusion will be that it is the growth 
in norm of the coefficient vector that in practice limits the achievable error, so this norm should be kept as 
small as possible to retain high accuracy. The reader should take care throughout not to confuse statements 
about the coefficient norm (which depending on the choice of MFS charge points may either grow or not 
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Fig. 1. Geometry for the MFS in the unit disc. 

grow with N as the error converges to zero) , with statements about the condition number of the problem 
(which always grows with N since the MFS (2) approximates a single-layer operator (3) which is compact). 

In Section 3.1 we move to general analytic domains, and review results for the analytic continuation of 
solutions u of (1), in particular how both the boundary data and the domain shape may lead to singularities 
in the continuation of u. In Section 3.2 we explore the use of the exterior conformal map in choosing the 
charge points for several more complicated and nonconvex domains. We propose and provide evidence 
for conjectures in general domains which arc analogous to the theorems on convergence rate and stability in 
the unit disc. In Section 3.3 we propose and demonstrate a method for choosing charge points well-adapted 
to the singularity locations and the wavenumber k, that outperforms the conformal mapping method by a 
large margin. 



2. The MFS on the unit disc 

In this section we analyse the accuracy and coefficient sizes that result in the unit disc = {x : |x| < 1}, 
for the MFS using charge points = Re l ^ j , j = 1, 2, . . . , N, with <j)j = 2irj/N, that is, equally spaced on a 
larger circle of radius R > 1. See Figure 1. We identify K 2 with C. 

Before we embark we need the definition of the Fourier series for a function g 6 L 2 ([0, 2tt]), 

g(9) = J2 9(™y m<> , 9(m) = 5- / 9{&)e~ ime dQ. (6) 

— rv\ w 



771 — ~ OO 



Parseval's identity is then 

oo 

\\9\\U[o,2n]) = 27r E \9{m)\ 2 =:27r\\g\\% (z) . (7) 



m— — oo 



We also need to represent the coefficient vector a := {ctj}j=i....,N in a discrete Fourier basis labeled by 
-N/2 < k < N/2 (we will always choose N even), 



N/2 1 AT 



E «* e "*'> a fe = ^|>e-^, (8) 

k=-N/2+l j=l 

where inversion follows from Y^,jLi e 2mk ^ N = N8^' with the periodized Kronecker delta defined by 
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a) wavenumber k=7, R=2 



b) wavenumber k=500, R=1.2 
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Fig. 2. Comparison for the unit disc of layer-potential eigenvalue magnitudes |s(m)| given in (12) against various asymptotic 
expressions: 'Laplace' (14), 'improved' (44), and 'uniform' (45). a) low wavenumber, b) high wavenumber. Panel c) shows 
density plot of matrix elements (20) of Q for N = 10, in the domain \m\ < 30. In c) we chose unrealistically small values of R 
and k in order to make the super- and sub-diagonals more visible. 



sir 



1, k = j (mod iV) 
0, otherwise. 

2 , . . 



(9) 



Parseval's identity now gives \a\ 2 — o: | 2 , where |a| := (|a?i 
norm. 



1/2 

|a„| 2 ) is the standard Euclidean 



2.1. Map from layer potential to Fourier basis on the unit circle 



For simplicity we first consider the layer potential version of this problem, which can be interpreted as 
the TV — > oo limit of the MFS. The single layer potential lying on the outer circle T = {y : \y\ = R} is 

«(x) = l - J 2 J H^(k\K- Bj+\)g{<f>)d4. (10) 

Note g € i 1 ([0, 2tt]) is the density with respect to angle measure d<j) rather than the usual length measure 
Rd<t>. The Fourier-Bessel decomposition of a fundamental solution located at Re 1 ^ ', evaluated at x = re is, 
using Graf's addition formula [1, Eq. 9.1.79], 

l -H^\k\^-Re^\) = \Y. HHHkR)J m (kr)cosm(6-<f>) = % - £ (kR)e~ m ^ ■ J m (kr)e lme (11) 

where the second step involved the reflection formulae [1, Eq. 9.1.5] J- m (z) = (— l) m J TO (z) and Y- m (z) — 
(— \) m Y m {z). Hence the Fourier-Bessel coefficients are jHm\kR)e^ tm ^ . The restriction of (10) to x e dfl 
gives a single-layer operator S : L 2 ([0,2tt}) — > L 2 ([0,27r]) of convolution type, which is therefore diagonal 
in the Fourier basis {e lm6l } me z and entirely described by its eigenvalues. Comparing (10), (11) and using 
orthogonality gives u(m) — s(m)g(m) where the eigenvalues of S are 

s(m) = ^H^(kR)J m (k). (12) 



Remark 1 Since the Hankel function (real argument) is never zero, an eigenvalue can vanish only when 
Jm(k) = 0, corresponding to a Dirichlet eigenvalue (resonance) of Q. In contrast if Y a were chosen as 
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the fundamental solution in (2), Y m (kR) may accidentally be very small giving poor or spurious numerical 
results, although in practice this happens rarely [7j. In general one may avoid this problem by using (i/4)(Y"o + 
inJo), for any real n ^ (an analogous idea is used in layer potentials, p. 48 of [4]). 

Since its kernel is continuous S is compact, so limimi^oo s(m) — 0. The compactness of S means its inverse 
is unbounded, and we expect to find arbitrarily large ||g|| needed to represent certain unit-norm boundary 
functions v. 

In the Laplace (k — > 0) limit we recover the following known result (e.g. [13, Eq 3.2]). We use the small- 
argument asymptotics J m (k) ~ {k/2) m /T(m + 1) and Y m {kR) ±T{m){kR/2)-' m for all integer m > 0, 

and the reflection formulae, and get 

s(m) -» ^|#" H 5 mGZ\{0}, fc^O. (13) 

To analyze convergence rate we will need the asymptotic behavior as \m\ — > oo for fixed k. Using in (12) the 
large-order asymptotics (9.3.1 in [1]) J m (z) <~ -^==(ez /2m) m and Y m (z) ~ —^J -^(ez/2m)~ m , where z is 
fixed, gives the leading-order behavior 

s(to) - ^-|i?- |m| , |m|-»oo, (14) 

which coincides with the Laplace case (13). We therefore have the exponential uniform bounds, for some 
constants c s and C s depending only on R and k, 

^i?-!" 1 ! < \s(m)\ < QjR-W < C s R- m , m e Z\{0}, (15) 
\m\ \m\ 

where C s is chosen large enough such that also |s(0)| < C s . Figure 2a shows that for low wavenumbers the 
leading-order asymptotic is reached rapidly, hence the smallest possible ratio C s /c s is not too large (here 
about 10 3 will suffice). This asymptotic becomes accurate well before the dynamic range begins to exceed 
machine precision (e mach w 10~ 16 for double precision). However, the situation can differ radically for large 
wavenumbers, as Figure 2b illustrates. Here the Laplace asymptotic is not relevant for the eigenvalues within 
a factor e mach of the maximum. Worse still, the ratio C s /c s must be exceedingly large (several tens of orders 
of magnitude). We will present more useful asymptotic approximations for the eigenvalues in Section 2.4. 
For now we need only (15) to prove exponential convergence rates. 

2.2. Map from MFS coefficients to Fourier basis 

We now adapt the above to the discrete source case. Define the density 

AT 
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(# = 5>^-&). (16) 



It follows that 



u W (e *) = if^ajH^ik^ -Re^\) = [Sg){cj>). 



We have 



\\Sg\\h ao ,2n])=' 2 M\S9\\(2 { z)= 2 ^ \K'm)g{m)\ 2 ^ — ^ \s(m)a m mod w ' 2 



2tt 

m— — oo 



where m mod TV denotes the unique integer lying in the range —N/2 + 1, . . . , N/2 which differs from m by 
an integer multiple of N. The last equality follows from the Fourier series representation of (16), 

1 N N 
= — J^aje = —a m modN, m £ Z. (17) 

J = l 
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Applying Holder's inequality and (15) we obtain 

N 2 °° N 2 R 2 4- 1 

^ 2 2 m— — 00 

Define the operator Q : l w — > £ 2 (Z) by := Sg. Therefore Q maps the discrete Fourier coefficient vector 
a. to the Fourier series coefficients on the boundary Oft. From (18) we immediately obtain the following. 
Lemma 2 For R > 1 the operator Q is bounded. Furthermore, 



N R 2 + l 

w/iere ||Q|| := max d£K iv\ {0} ll9 "^ 2 ( z > ■ 

The action of Q is that of a generalized matrix of width TV but (bi-)infinite height, 

N/2 

u(m) = ^ Qmk&k, for m e Z. (19) 

fe=-JV/2+l 

From Qat(m) = (Sg)(m) = s{m)g{m) and (17) it follows that the matrix elements are 

qmk = ^8{m)5^. (20) 

Fig. 2c) shows a greyscale picture of a piece of the resulting matrix Q. Notice that it is dominated by a main 
diagonal proportional to the diagonal of the S operator defined in (12), but with (exponentially) smaller 
entries on an infinite sequence of super- and sub-diagonals. This off-diagonal part can be interpreted as 
aliasing 'overtones' due to discrete sampling of a continuous layer potential. In the Laplace case using (13) 
in (20) recovers the results of Katsurada [15, Lemma 1, case 2]. 

2.3. Convergence rate and coefficient sizes in the disc with analytic data 

We are now in a position to express the boundary error norm (4) in terms of a. Combining with (7) and 
(19) gives 

t[d] ^V2^\\Qa-vy m , (21) 
where v e £ 2 (Z) is the p a recipe for locating charge points which allows us to reach error norms of typically 
10~ n on a wide variety of analytic domains. At high frequencies of order only 3 points per wavelength are 
needed, which compares very favorably to boundary integral methods. 

vector of Fourier coefficients of the boundary data v on the unit circle. Assume that v can be analytically 
continued to the annulus {z e C : j < \z\ < p} for some p > 1, that is the closest singularity of the analytic 
continuation of p has the radius p or 1/p. We then have asymptotically exponential decay of the Fourier 
coefficients, 

\v(m)\ - Cp Hm| , \m\ 00, (22) 
for some constant C. A simple example is boundary data arising from an n th -order pole v(z) = Re (z — p)~ n 
for z € <9f2, n = 1, 2, . . . . 

Minimizing (21) over d is a least-squares problem involving the generalized matrix Q. But since the 
columns of Q arc orthogonal this separates into N independent single-variable minimizations. We may use 
a diagonal approximation to choose a which is sufficient for the following convergence rate bounds. 
Theorem 3 Let R > 1 and N be even. For analytic boundary data v obeying (22), the minimum boundary 
error (4) achievable with the MFS in the unit disc satisfies 



t < < 



Cp- N '\ p < R 2 

C^NR- N , p= R 2 (23) 
CR~ N , P >R 2 
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where each time C means a different constant which may depend on k, R, and v, but not N. Furthermore 
if v is analytically continuable to an entire function, the last of the three cases holds for any R > 1 . 
Proof: We choose coefficients a m — v(m)/q mm for —N/2 < m < N/2. This exactly matches the Fourier 
coefficients in this interval, therefore errors are due only to frequencies lying outside the interval. (21), (20) 
and the triangle inequality in £ 2 (Z) give 



1/2 



where 



and 



t = 



El = 



2ir ]T \{Q*){m) - v(m)\ 2 
m £[-f+i,f] 



< V2n{E u + E v ), 



E kq«)mi 2 

m£[-# + l,f] 



- E 



v(n) 



s(n) 



(24) 



El 



\v(m)Y 



(25) 



E 

m£[-#+l,§] 

We can bound both error terms since all terms in the sums have exponential bounds. First we note that 
using (15) and (22) gives 



<Ci 



-#<«<# 

E 

-f<«<f 

n^O 



2\n\ 



n 



b/0 



|&iV + n| 2 



fi(O) 



8(0) 



E 

f>#0 



|WV| 2 



2\n\ 



E^ 2|bJV+n| + c 2 

6#0 



t>(0) 



5(0) 



E^ 2|bJV| 

6^0 



(26) 



for sufficiently large constants C\ and C2. We can bound 



^ ~ 2 ~~ 2 

Mo 



(27) 



for a large enough constant C3. Inserting (27) into (26) and absorbing 

, 2|n| 

P 



fi(0) 



§(0) 



into the constants gives 



El < CR~ 2N 



(28) 



for a sufficiently large constant C. Similarly E 2 < Cp~ N follows from (22). We now study the sum in (28). 
For p < R 2 , (28) can be estimated by Cp~ N for some constant C > 0. This means both E 2 and E 2 have the 
same exponential decay Cp~ N . For p = R 2 , the sum contains N equal terms and therefore E 2 < CNRr 2N , 
which decays slower than the bound Cp~ N on .E 2 . For p > R 2 , (28) can be estimated by CR~ 2N for some 
C > 0, which means is of higher negative order in N than and can be dropped. For the case of v 
continuable to an entire function, we may take p — > 00 and the case p > R 2 applies. □ 



This is a generalization of a result of Katsurada [12] from the Laplace to the Hclmholtz problem. Since 
only the exponential bounds (15) and no other information about §(m) was used, the convergence rates are 
identical to those for Laplace with the same boundary data. 
Remark 4 An interpretation of the two main convergence rate regimes is, 

- v is 'not relatively smooth' (R > ^fp, i.e. 'distant' charge points): errors are limited by the absence of 
Fourier modes beyond a frequency N/2 in the MFS basis, hence rate is controlled by the boundary data 
singularity p. 

- v is 'relatively smooth' (R < ^fp, i.e. 'close' charge points): errors are limited by aliasing errors due to 
the discrete representation of the single layer potential, hence rate is controlled by R. 
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Remark 5 By keeping track of the constants in the above proof one can check that C is at least as big as 
C s /c s , which, as discussed, must be very large for large k. Thus while the convergence rates in the theorem 
must be reached asymptotically as N — > oo in exact arithmetic, we cannot expect the bounds to be numerically 
useful in practice at large wavenumbers. 

We are interested in how the coefficient norm \a\ grows as we reduce the boundary error in the MFS 
representation (2) by increasing N. Firstly, it is easy to show that when the MFS charge points are closer 
than the nearest singularity, the coefficients need not grow. 

Theorem 6 Let be the unit disc, and R < p, with fixed analytic boundary data v obeying (22). Then as 
N — > oo there exists a sequence of coefficient vectors a with bounded norm \a\, with corresponding boundary 
error norm (21) converging as in Thm. 3. 

Proof: We choose coefficients as in the proof of Thm. 3, which therefore give the desired convergence rate. 
Then using (22), (20) and (14), 

vim) 2irv{m) Jml / i?\ M C ( i?\ M N N 



q mm N s(m) N\pJ -2\pJ ' 2 "2 

for some constant C. For R < p this is an exponentially decaying sequence so, independent of N, \a\ is 
bounded by a constant. □ 



More problematically, the coefficient choice used in the above two proofs would then imply in the case 
p < R that \a\ diverges exponentially with TV. However, it is not immediately obvious whether there is a 
different choice of a m which avoids exponential growth. The following theorem excludes this possibility by 
showing that when the singularity in the analytic continuation of the boundary data is closer than the MFS 
source points, any convergent sequence of coefficient vectors a must diverge in norm in this way. 
Theorem 7 Let f2 be the unit disc, with R> p. Let the boundary data Fourier coefficients decay no faster 
than (22), that is, for some constant c v , 

\v(m)\ > c v p-\ m \. (30) 



For any positive even N satisfying N > 3 



N min , where N„. 



2 max 



In 



it V** 1 



let a. be 



coefficient vector such that the MFS representation (2) has a boundary error norm (21) satisfying 



t < cp 



-N/2 



where c is a constant independent of N . Then 



N/2 



(31) 



(32) 



for some constant C which may depend on k, R, and v, but not N. 

Note that (31) is the appropriate convergence rate for the case p < R 2 derived in Theorem 3. 
Proof: For even N fix N > iV min + 3. Using (20) in (21) implies the trivial bound 

N 



2tt 



s(m)a m mod at - v(m) 



<t, 



for all m € Z. 



(33) 



Define the (positive) maximum Fourier frequency F := ^ — K, where K is the unique integer such that 



N„ 



< K < 



1. 



2 - 2 

Note that K is independent of N. One can verify using (31), (30) and the definition of 7V min that 



v(m)\>\ -t, for all \m\ < F. 



(34) 



In the frequency range \m\ < F, it follows from (33) and (34) that v(m) is sufficiently large relative to t to 
bound the coefficients away from zero, 
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Fig. 3. Convergence and coefficient sizes as a function of N, for the MFS approximation to the interior Helmholtz BVP in 
the unit disc given boundary data corresponding to a single source (39) at radius p. The wavenumber is low (k = 8). The 
MFS sources are at R = 1.5. For visual comparison the relevant power laws from Theorems 3 and 7 are shown (sometimes the 
constants have been chosen to match the data). There were M = 240 boundary points. 



> 



2-7T — t/y/2n n \v(m)\ 



> 



N \s(m)\ ~ N \s(m)\ ~ C S N 

where the last step used (15) and (30). Choosing the maximal frequency m 



a F 



> 



TTCv 

c7 



K 

N 



> 



TTC V 

2Cl 



< \m\ < F (35) 
F = ^ — K we obtain 

' (36) 



^Vmin + 3 



Here the latter inequality follows from 

1 K 

2 ~ N 



> 



N min /2 + 1 



N min + 3 2(N min + 3) 

Absorbing the ./V-independent factors of (36) into a constant and noticing that the Euclidean norm of a 
vector is at least as large as its largest component we have 



— — — 'R\ N/2 

Vn\ol\ > Vn\& f \ > cVn[ — \ 



■ 



for a sufficiently small constant C > 0. 



(37) 

□ 



An immediate consequence is that if p is known for given boundary data, to prevent exponential growth 
in coefficients one should restrict the charge point radius to R < p. To illustrate this and Thm 3, we finish 
this section with some numerical experiments at low wavenumber. The implementation was standard, as 
follows. The integral in (4) is approximated using uniform quadrature with M equally-spaced boundary 
points {x m } m= i...M- Specifically the M-hy-N matrix A has elements 

A mj := l -H^\k\^- yj \), (38) 

and the boundary- value vector v £ C M has elements v m :— n(x m ). The resulting linear system (usually 
overdetermined, M > N, in our work) Aa = v was solved in the least-squares sense via the QR decompo- 
sition (MATLAB's backslash command) in double-precision arithmetic. The boundary error norm then is 
approximately t — y/\d£l\/M \Aa — v\. 
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Fig. 4. MFS approximation to the interior Hclmholtz BVP in the disc with given boundary data corresponding to a single 
source (39) outside the domain at p = 1.2, with TV = 80 and R = 1.4. a) Re u^ N \ b) residual Re — u) on a color scale 

3 X 10 4 times more sensitive than a). 

In Fig. 3 we show convergence of t using boundary data 

v(z) = --Y (k\z- p\), zedn (39) 

with real p > 1, that is, a single real- valued fundamental solution. The three panels illustrate the three 
cases of Theorem 3. In a) p > R 2 thus convergence is determined by R. In b) R 2 > p > R so we have 
transitioned to a convergence rate given by p. In both these cases the coefficient size \a\ is very close to 
constant (note by contrast that the condition number of A is growing) . However in c) p < R so convergence 
rate is again determined by p, but now \a\ grows exponentially at precisely the rate indicated by Theorem 7. 
Fig. 4 shows the resulting approximate field U W in the case c) , and the error function Notice that 

the error is oscillatory at Fourier frequencies of about N/2 (as predicted by Remark 4; this can be seen by 
comparing the alternating signs in Fig. 4b) to the angular spacing of source points) , is concentrated on the 
side of dfl nearest the singularity, and decays exponentially inside the domain (it is evanescent). We have 
also substituted v(z) = Re(z — p)^ 1 and find the convergence rates in Fig. 3 are very similar. We note that 
in each plot in this figure, the convergence eventually stops, as we now explain. 

2.4. Minimum achievable error in the disc for low and high wavenumbers 

So far we have proven results which hold in exact arithmetic. However machine precision limits the dynamic 
range of eigenvalues that may be used: since the MFS trial functions in (2) have typical size of 0(1) in fl (for 
any reasonable wavenumber), each coefficient otj will result in round-off errors of size roughly e mach aj in the 
numerical approximation u^ N \ Thus we expect convergence to stop when t reaches of order e mach times the 
coefficient norm \ac\. This behavior is well illustrated in the three plots of Fig. 3: convergence stops when the 
ratio between t and |a| reaches roughly 10 -16 . In c) the coefficient growth thus limits achievable error norm 
to only about 10 -5 . Such premature halting of convergence has been observed in the Laplace (k = 0) case 
in the disc [12,16,28] but not analyzed much before. We analyse this in the Hclmholtz case after observing 
the following consequence of Thm. 3 and 6. 

Remark 8 For any wavenumber, for boundary data with a given singularity radius p, the choice of charge 
point radius R in the range ^fp < R < p leads to both optimal asymptotic convergence rate t ~ p~ N I 2 and a 
lack of coefficient growth. 

It is useful to have a heuristic model which predicts, for general p and R in the unit disc, both the 
lowest achievable error norm and the basis size N required to achieve it. We spend the rest of this section 
constructing then testing such a model. We first consider low wavenumbers, that is, ones where the Laplace 
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Fig. 5. Convergence and coefficient sizes as a function of N, for MFS approximation to the interior Helmholtz BVP in the disc 
at high wavenumber k = 500. The boundary data corresponds to a single source (39) outside the unit disc at radius p, with 
MFS sources at R. For a) p = f.f, R = 1.2, b) p = 1.01, R = 1.05. The 'predicted' curves are given by (42) using (41), with 
(45) modeling both s(m) and v(m). 




Fig. 6. a) (left plot) Minimum achievable boundary error norms to and b) (right plot) corresponding basis sizes No, predicted 
for a selection of source radii p and MFS source point radius R, for the Helmholtz BVP in the unit disc at high wavenumber 
k = 500. The predictions are done in the range R > p using the model in Sec. 2.4. Note that in b) all the graphs for different 
p lie on top of one another. 

asymptotic form (13) is relevant for the relevant eigenvalues (those no smaller than e mach times the largest 
eigenvalue). We discuss both the case R < p where no coefficient growth occurs, and the case R> p when the 
coefficient norm grows. The proof of Theorem 7 suggests that, at least when s(m) is exponentially decaying, 
the diagonal approximation for the MFS Fourier coefficients (see (29)) 

2n vim) 
JS s(m) 
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approximates the true least-squares Fourier coefficients well apart from an O(l) number of them lying at 
extreme frequencies near to = ±N/2. For R > p these Fourier coefficients grow exponentially, so dropping 
an 0(1) factor we may consider only the largest coefficient's contribution to the ^ 2 -norm, namely that at 
to = N/2. For R < p there is no growth hence the norm is dominated by the coefficients of size O(l) at 
to s=a 0. Therefore we have the order-of-magnitude estimate at a given N, 

where the maximum- value operation combines the two cases. Similarly, since the boundary data coefficients 
die exponentially, following Remark 4 and ignoring an 0(1) factor we may suppose 

t*\v(N/2)\. (42) 

We define No to be the N at which convergence stops, for this we use the round-off error consideration 
t/\ct\ ~ e mach discussed above. In the case R < p this implies that convergence halts when t reaches of order 
Emach! as observed in Fig. 3a, b. However for R > p, combining (41) and (42) we get an implicit equation for 

\/N~o\s(No/2)\ e mach (criterion for halting of convergence, for R > p) . (43) 

The minimum achievable boundary error is then given by (42) with the substitution N — No- As an 
illustration, using the (Laplace) asymptotic form (14) for the eigenvalues approximately predicts (dropping 
algebraic factors) that iVo ~ 2 ln(l/e mach )/ lni?. For the parameters of Fig. 3c this gives No s=a 180, then 
using (22) with C = 1 gives to ~ 1CP 4 which, given the heuristic nature of our model, agree well with the 
observed behavior. 

We now briefly discuss the case of high wavenumbcr. In Fig. 2 we saw that the (Laplace) asymptotic form 
(14) is not useful for predicting relevant eigenvalues at high k. We may derive (Appendix A) the asymptotic 

S(m) ~ J_ R ^n\ e k 2 (B 2 -l)/4m , , ^ (44) 

v ' 2\m\ > i i j \ ) 

which Fig. 2a), b) shows is a much improved approximation, but still not useful for the relevant eigenvalues 
at k = 500 (or beyond). Therefore we use the WKBJ method to derive (see Appendix A) a uniform 
approximation for the eigenvalue magnitudes, defining a 2 = m 2 — j, 



\s(m) 



[(k 2 -a 2 )(k 2 R 2 -a 2 )f V4 , m<k 

]-[(a 2 -k 2 )(k 2 R 2 -a 2 )Y 1/A e^ k \ k < m < kR (45) 

I [(a 2 - k 2 )(a 2 - fc 2 i? 2 )]- 1/4 e J »( fe )- J »( feii ), m>kR 



where 



I a (x) := Va 2 -x 2 -a\n[(a+Va 2 -x 2 )/x}. (46) 
More precisely this is an estimate of the amplitude in oscillatory region (to < k) of J m , and the absolute 
value in the evanescent region (note can also be complex oscillatory but its magnitude never is). In 
the oscillatory region individual s(to) values cannot be predicted: rather they are distributed in the range 
[—1,1] times the approximate amplitude (45). Fig. 2 shows this is a highly accurate asymptotic form in all 
regions apart from the two turning-points (J m (fc) is at its turning point for to th k whereas for Hm\kR) 
this occurs at m w kR). The estimate has algebraic singularities at these two turning-points, but they are 
weak enough that it is still useful. 

Finally, we compare numerical convergence results against this model at high wavenumber. Fig. 5 shows 
convergence and coefficient norm at k = 500 (about 170 wavelengths across the domain) for boundary data 
(39) deriving from an exterior fundamental solution. Its boundary data Fourier coefficients v(m) are given 
by the same formula (12) as the MFS eigenvalues (and hence the same approximation (45)) but with the 
substitution p for R. To compute the curves shown as 'predicted', we used this approximation in (41), and 
(42) to predict the error norm. It is clear that, up to the point when convergence halts, the predictions for 
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both error norm and coefficient norm are very close to observations (the largest deviations being spikes due 
to algebraic singularities discussed above; in a) these are at N = 1000 and 1100). 

A crucial common feature is that no convergence happens until N — 2k, since s(m) remains large for 
\m\ < k (see Fig. 2b). One interpretation of this is that 2k, corresponding to 2 degrees of freedom per 
wavelength on the perimeter, is the Nyquist sampling frequency for /c-bandlimitcd functions on <9f2; in 
physics this is known as the the 'semiclassical basis size' [2]. In panel a) of the figure, p = 1.1 so the 
singularity is k(p — l)/2ir ss 8 wavelengths from the boundary. In this case convergence is rapid, dropping 
ten orders of magnitude between N = 1000 and N = 1150. Convergence then halts (compare (43) which 
predicts Nq 1175 and to ~ 3 x 10~ 14 ). The number of boundary quadrature points was M = 1500 in a) 
(only 3 points per wavelength). This can be chosen to be so small since boundary functions v and u^lan 
have exponentially-decaying Fourier coefficients beyond frequency 2k, giving spectral convergence. In panel 
b) p = 1.01 and R = 1.05, giving both slower convergence and growth in coefficient norm. The predictions 
Nq Ri 1708 and to s=a 8 x 10~ 7 are again reasonably close to observations. 

How can R best be chosen to achieve the lowest boundary error for a given high wavenumber k, and pi 
We use the above model to compute Nq and hence to for a variety of p and R at k — 500, in Fig. 6. Here the 
smallest p — 1.002 corresponds to a singularity 0.16 wavelengths from the boundary. The conclusion is that 

~ £mach appears to be always achievable as R tends to p from above, as expected from Remark 8; however, 
the basis size required to do this diverges as p — > 1 + . b) also shows that there is a limiting basis size of 
about N ;=y 1180 (not much larger than 2k) for which arbitrarily large R may be used, but with this choice 
to becomes 0(1), hence not useful, as p — > 1 + . In conclusion, we may state that in the high-wavenumbcr 
limit, if the nearest singularity in the boundary data is at least a few wavelengths away, then both the basis 
size N and the number of quadrature points M can approach 2 per wavelength while achieving an error 
close to machine precision. 

3. The MFS on analytic domains 

In this section we present results for the MFS on arbitrary analytic domains. On the circle we have shown 
that the MFS coefficients start growing exponentially if the radius R of the charge points becomes larger 
than the distance p of the singularity. In this section we demonstrate that also on general analytic domains 
the position of the charge points relative to the singularities of the analytic continuation is crucial for the 
accuracy and numerical stability of the MFS. 

3.1. Analytic continuation of solutions 

The question of analytic continuation is to find a domain fi D D, and a function u such that Am + Am = 
in f2 and u\q — u. Since solutions of the Helmholtz equation are real analytic it follows immediately that u 
is unique. 

A classical result of analytic continuation is reflection on a straight arc T, on which u satisfies tt|r = 0. 
Without restriction let T be a subset of {iy : y € M}. Then u can be continued across Y by setting 
u(—x,y) := —u(x,y) (see also [5]). In [9] Garabedian extended these results to the case that Y is an 
arbitrary analytic arc for which Mp = 0. More general reflection principles for linear elliptic PDEs of the 
type Aw + a{x,y)u x + b(x,y)u y + c(x,y) = 0, where a(x,y), b(x,y) and c(x,y) are real analytic functions 
were treated by Lewy in [23]. He stated his results for arbitrary Dirichlet, Neumann and mixed boundary 
conditions but restricted T to be a straight line. Representations of the analytic continuation for the case 
that T is not a straight line were given by Millar in [24]. In [25] he discussed more in detail the analytic 
continuation of solutions of the Helmholtz equation. 

Millar shows that there are two possible sources for singularities of the analytic continuation u of u. The 
first one comes from singularities of the analatic continuation of the boundary data /. The second possible 
source of singularities is introduced by the shape of <9£1 Let Z(s) = x(s) + iy(s) be a parameterization of 
dCl, where s £ [0, 2ir]. Assume that x(s) and y(s) are real analytic and that |Z'(s)| ^ in [0, 2n}. Then there 
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Fig. 7. Illustration of map Z(s) and its inverse S(z) defining a boundary curve dQ. The Schwarz function involves composition 
of Z(s) and S(z). 

exists a complex neighborhood of [0, 2tt], in which Z{s) is holomorphic and invertible. We denote its inverse 
by S(z) and define the Schwarz function 

G(z) := Z{S{z))=Y{S{z)). 

Millar showed that except for special cases the singularities of G(z) outside f2 are also singularities of the 
analytic continuation u of u. 

The Schwarz function has been studied in [6] . It is independent of the parameterization of dQ and has an 
interpretation in terms of reflection principles on analytic arcs. Assume that z\ is a point close to dfl. Then 
its reflection on dil can be obtained by the following steps (see Figure 7). 

(i) Compute t\ = S{z\). 

(ii) Reflect t\ on the real line to obtain the point t 2 := t\. 

(iii) The reflection Z2 of z\ at dQ is now obtained as 

Z2 = Z(t 2 ) = Z{tl) = Z(S(zij) = Gjzij. (47) 

Fig. 8 shows the singularities of the Schwarz function on three different domains, a rounded triangle, an 
inverted ellipse and a crescent. The domains are defined with default values for the parameters a\ through 
a4, as follows. 

Rounded triangle: Zt{s) = e ls + aie _2,s , a± = 0.3 
Inverted ellipse: Zie(s) = 1+ „ , a-i = 0.25 

Crescent: Z c {s) = e is - a 3 = 0.1, a 4 = 0.9 

Branch type singularities are denoted by '+' and pole type singularities by '*' in Figure 8. The branch 
singularities in all three domains are of square root type (see [25] for an analysis of the branch behavior of 
G(z)). The crescent has exterior singularity of pole type at z — — l/aj. For the interior Hclmholtz problem 
only the exterior singularities of G are important since these are points where, for generic boundary data, 
the analytic continuation u of u becomes singular. Conversely if we had an exterior Hclmholtz problem then 
the interior singularities would determine the singularities of the analytic continuation. 

3.2. Using exterior conformal map to place the charge points 

A natural generalization of the MFS on the unit disk to general analytic domains can be defined in terms 
of the conformal map from the exterior of the unit disk to the exterior of the domain. This was investigated 
in the Laplace BVP case by Katsurada [14]. 

Let n be a simply connected domain with analytic boundary dfl. We can parameterize <9f2 using the 
exterior conformal mapping function 
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Fig. 8. Domains a) rounded triangle, b) inverted ellipse, c) crescent. Branch singularities of the Schwarz function are denoted 
by '+' and pole type singularitcs by '*'. Note for the crescent, the '+' signs are inside f2 but very close to the boundary. 

z = V(w) = cw + c + — + —Z- + . . . , c>0 

w w l 

which maps the exterior D\ :— {w : \w\ > 1} of the unit disk to the exterior of 0. The quantity c is called 
the capacity of ft. We denote the inverse map by w = $(z). Since dfl is analytic ty(w) can be analytically 
continued to a domain D r :— {w : \w\ > r} for some < r < 1. We denote the conformal radius of a point 
z € C\Q by p x := \&(z)\. 

In the notation of the previous section we may write this parametrization as Z(s) = ^(e ls ) since the 
unit disc is parametrized by w — e ts . Using that the reflection of a point z on the unit circle is given 

by z' — = the Schwarz function may now be written G(z) = *(l/<f>(z)); it follows that it is analytic in 
{zeC: 1 < Pz < i}. 

In the unit disk case we placed the MFS points equally distributed on a curve with radius R. For general 
analytic domains we now place the points on a curve := {z : p z — R] with constant conformal radius 
p z = R. On this curve we distribute the points equally spaced in conformal angle, that is 

y . := ^(e^/JV), j = l,...,N (48) 

If fi is the unit disk this definition coincides with that of Sec. 2. Replacing the disk radii R and p in Theorem 
3 by the corresponding conformal radii we obtain the following conjecture for the rate of convergence of the 
MFS for Helmholtz problems on general analytic domains. 

Conjecture 9 Let t be the error of the MFS as defined in (4) by placing the MFS points equally distributed 
in conformal angle at a conformal distance R around f2. Let p > 1 be the conformal radius of the closest (in 
the sense of conformal radius) singularity of the analytic continuation of u. Then 

f Cp- N l\ p < R 2 , 
t < I (49) 

" \CR- N , P >R 2 , 

where C is a constant that may depend on fi, k, R and v, but not N. Furthermore, if u continues to an 
entire function, the latter case holds for any R > 1. 

Remark 10 The first case of this conjecture was proved in the Laplace case by Katsurada in [14] under 
additional restrictions on the analytic continuation of ^ into the unit disk. In numerical studies we have 
observed that these conditions are not necessary to achieve the given convergence rates, so do not include 
them in our conjecture (compare also Remark 3.2 of [14]). We do not state a conjecture for the case p = R 2 
since numerically it cannot be established if for general domains the same algebraic factor is needed as for 
the disk. 

In Figure 9 we plot the observed error t for the MFS on the inverted ellipse of Figure 8, for wavenumber 
k = 5 and constant boundary condition v = 1. The estimated rates from Conjecture 9 are denoted by dashed 
lines. The three plots correspond to MFS points placed at the conformal distances R = 1.03, R = 1.12 and 
R = 1.2. The conformal radius R = 1.12 is also the approximate conformal radius p of the two singularities. 
The corresponding MFS curves are shown on the right of Fig. 9. The estimated convergence rates are in all 
three cases in good agreement with the observed error t. 
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R=1.03 R=1.12 R=1.2 




Fig. 9. Left three plots: Estimated (dashed lines) and observed ('+') rates of convergence of the MFS on the inverted ellipse for 
charge points with conformal radii R = 1.03, R = 1.12 and R = 1.2. The corresponding coefficient norms |oe| are denoted by 'o'. 
Right plot: Domain and positions of the charge points (shown as dots) for the above three conformal radii, with branch-type 
singularity locations ('+' signs). 

In the disk case we can observe exponential coefficient growth once the radius R of the charge points is 
larger than the radius p of the singularity of the analytic continuation of u (see Thm. 7) . For general analytic 
domains we observe a similar behavior. To demonstrate this we plotted in Figure 9 also the norm \ct\ of the 
MFS coefficients. As long as the conformal radius R of the MFS points is smaller or equal to the conformal 
radius p of the singularities we do not observe any growth of jo:] for growing N (first two plots of Fig. 9). In 
the third plot we have R > p and |a| grows exponentially for growing N . It is instructive to compare this 
figure panel by panel against Fig. 3. 

This leads to the following conjecture, which mirrors Theorem 6 and 7. 
Conjecture 11 Let the MFS charge points be chosen equally spaced in conformal angle on T r := {z : p z — 
R}, the curve of all points with given conformal radius R > 1. Let p be the conformal distance of the closest 
singularity of the analytic continuation of u and let a be the vector of coefficients of the MFS basis functions 
that minimizes t. We have 

\a\ > C-i N 

for some 7 > 1, and a constant C that does not depend on N, if and only if R > p. 

In numerical experiments we observed also for other types of MFS curves that there is only coefficient 
growth if the curve encloses a singularity of the analytic continuation. Hence, a more general conjecture can 
be stated (similar to results known for the scattering case [21]). 

Conjecture 12 Let T be any Jordan curve enclosing fl, with dist{T, dfl) > 0, on which MFS charge points 
are chosen asymptotically densely. Then the coefficient norm \a\ that minimizes t grows asymptotically 
exponentially as N — > 00 if and only ifT encloses a singularity of the analytic continuation of u. 

In Figure 10 we show the coefficient norm |a| and the approximation error t for a growing conformal 
distance R of the MFS source points and fixed number N in four different cases: the unit disk and the three 
domains from Figure 8. In all cases we have used k = 5. For the disk the boundary data is given by (39) with 
singularity location p = 1.2, and in the other three cases by v(z) = 1 (recall that here the Schwarz function 
introduces singularities in u). The vertical solid lines denote the conformal radius p of the singularities of 
the analytic continuation of u and the vertical dashed lines denote the square root p 1 / 2 . Since the Schwarz 
function for the rounded triangle does not have any singularities in the exterior of the domain the solution 
u can be analytically continued to an entire function. 

For the disk, the inverted ellipse and the crescent the error t does not decrease further once R passes the 
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Fig. 10. Growth of coefficient norm |ct| (circles), and least-squares approximation error t ('+' signs), for a) the unit disc with 
boundary data (39), and b-d) the three other shapes of Fig. 8 with constant boundary data v = 1. MFS points are chosen on 
curves with constant conformal radius R. The insets show domains, exterior singularities using same notation as Fig. 8, and 
curves with conformal radius R, = 1.24 and 1.49. The parameters are k = 5. N = 200 for a-b, TV = 400 for c-d. The number of 
quadrature points M on dQ is chosen sufficiently large throughout. 

dashed line. This can be expected from Conjecture 9 since the upper bound on the error t does not decrease 
any more for fixed N and R > p 1 ! 2 . 

For these three domains we can also observe exponential coefficient growth of |a| for fixed N when 
R > p. For the crescent this exponential growth already starts earlier. However, this is not a contradition to 
Conjecture 11. The conjecture treats the case of fixed R and N — > oo. This does not exclude the existence 
of transient growth effects for R < p. An explanation for these transient effects in the crescent case is that 
close to the pole-like singularity of the Schwarz function we need a very high number N of basis functions 
to sufficiently resolve a highly-oscillatory Helmholtz field. 

Another interesting special case is the rounded triangle. Since the analytic continuation of u is an entire 
function, by Conjecture 11 we do not expect any exponential growth of \a\. Indeed, Fig. 10b shows that \a\ 
stays virtually constant as R increases. 

In this section we have demonstrated with several numerical experiments that the behavior of the MFS 
for Helmholtz problems on general analytic domains is similar to the unit disk case if we choose the MFS 
points on curves with constant conformal radius R. We stated two conjectures about the approximation 
error t and the coefficient norm |a| as N — > oo. From our numerical experiments and under the conditions 
that the conjectures hold it follows that the optimal conformal radius R for the position of the singularities 
is given by R — p 1 ^ 2 , where p is the conformal radius of the closest singularity. This ensures a maximum 
rate of convergence for t while keeping \a\ from growing exponentially (as in Remark 8). Furthermore, in 
view of the results in the crescent case of Figure 10 it seems advisable to choose R not too close to p if 
the singularity is determined by a pole in the Schwarz function in order to avoid large transient coefficient 
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Fig. 11. Convergence rate for the crescent domain (with = 0.1, = 0.9) illustrating the concavity problem. The wavenumber 
is k = 3 and constant boundary data v = 1. For a) and b) MFS charge points are placed on an exterior conformal curve at 
R = yj~p where p = l/a4 is the singularity conformal distance. The point spacing is equal in a) conformal angle, and b) 
arclength. c) Adapted curve and spacing given by (50), see Sec 3.3. The dashed line shows the predicted convergence rate for 
a), <. 



growth. 



3.3. Using a singularity-adapted curve to place the charge points 

The above exterior conformal method has the following problem: in any concave parts of fl the exterior 
conformal map ^ has a very (in fact, exponentially) large gradient. This well-known property of conformal 
maps is related to the so-called crowding problem. This has two consequences for concave regions: the spacing 
of charge points according to (48) becomes very large, and as R is increased from 1, the curve Tr moves 
away from dfl very rapidly. Both these effects are illustrated by the MFS charge curve for the crescent in 
Fig. 11a. 

This means that Schwarz function singularities which are a moderate distance from a concave part of 
dft may actually have a conformal radius extremely close to 1. The net result is that if the coefficient 
growth of Conj. 11 is to be prevented, R must be very close to 1, hence by Conj. 9 the convergence rate 
is necessarily very poor. For example, in Fig. 11, the conformal radius of the pole in the Schwarz function 
is only p — \/a^ w 1.11, and the observed rate for case a) approaches the predicted p~ N / 2 (dashed line in 
the figure). One can attempt to fix the problem of the large point spacing by retaining the same MFS curve 
Tfl but choosing charge points equally spaced in arc- length, as illustrated in Fig. lib. Despite an initial 
improvement for small N, the asymptotic convergence rate turns out to be no better than in case a), and 
is believed to be the same (for errors t < 10~ 6 it performs worse, we believe due to a lack of point density 
near the 'spiked' parts of the crescent). 

In Fig. 11c we show a different 'adaptive' choice of MFS curve and point spacing which clusters charge 
points near the singularity but spreads them out (while taking them further away from dfl) away from the 
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Fig. 12. Generalized crescent analytic domain parametrized by (52), showing performance of adaptive charge points of Sec. 3.3. 
The top row is at k = 3, the bottom row at k = 100. a) and d) show the curve y(x) an d the locations (Xj > 3/(Xj )) m the 
.s-plane, the singularities s CT ('*' symbols), and the distance-limiting function |Z'(x)|/D max (dotted line), b) and e) show the 
charge locations (for clarity, TV = 90 has been used in both cases), c) and f) shows error norm convergence t ('+' symbols) and 
coefficient norm |cx| ('o' symbols). 



singularity. This gives a convergence rate over 5 times faster than the exterior conformal curve, for instance 
an error of t w 10~ 13 is reached for N = 140 as opposed to N — 730 for case a). This simple adaptive 
method is based on the idea of replacing the exterior conformal map by an annular one, as follows. 

From the discussion in Sec. 3.1, the map Z(s) defines an annular conformal map which is one-to-one for s 
in some strip around [0, 2ir], as in Fig. 7. The external singularities control convergence rate; we label them 
by a = 1, . . . , P, where P is the number of singularities. They have s-plane locations s a = Xa ~ i T a with 
T a > 0, Their minimum distance to the real axis is r := min CT r a . Katsurada et al. [16] have discussed using 
such an annular map to place charge points for the Laplace BVP, according to := Z(2irj/N — i log R) for 
some R > 1; note this is the annular map equivalent of (48). A related annular map method has been tried 
in the scattering case [11]. According to Conj. 12 in order to prevent coefficient growth in this case one would 
need to choose log R < r. This may be a severe restriction: for example one may check that the crescent 
formula Zq(s) in Table 3.1 is identical to this domain's exterior conformal map 'J, thus the concavity effect 
causes r to be very small. 

However, there are an infinite family of annular maps Z(s) which parametrize the same boundary c?£l, 
and these may have differing Schwarz singularity locations in the s-plane. Given an analytic domain Q one is 
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Fig. 13. Interior Helmholtz BVP solution u for the same domain as Fig. 12, at k = 100, with boundary data v = 1, using as 
basis size of TV = 525, and M = 1000 boundary points. Error norm is t = 5 X 10~ n . The boundary (thin solid line) and charge 
points ('+' symbols) are shown, and the field u from (2) is shown both outside and inside f2. CPU time was 2.9 s to compute 
the coefficient vector ct, and 5.4 mins to evaluate the solution u shown (2.4 X 10 points on a grid of spacing 0.005). 

free to choose between such parametrizations. Ideally, our goal is to choose one with as large a r as possible, 
to achieve a high convergence rate. However, we find it convenient to retain the given parametrization Z(s), 
and instead build a charge curve in the s-plane which no longer has constant imaginary part, and which 
captures the spirit of such a reparametrization. Our curve is given by s — x — iy{x), where x € [0, 27r] is the 
real part of s, and the function y is given by 

(50) 

where parameter values performing well in most domains are j3 = 0.7 (interpreted as a curvature factor), 
7 = 0.4 (expressing the curve's fractional distance to each singularity), and 

D max - max[l, 25/fc] (51) 

is interpreted as the maximum allowable distance of the curve from <9f2. Roughly speaking, (50) has the 
effect of bringing the curve close to d£l in the vicinity of each singularity (via each cos term in the sum), 
while allowing it to move up to -D max from the boundary in the absence of nearby singularities. Given the 
curve function y(x) a se t of N real values {xj} G [0, 27r] are then chosen such that their local spacing is 
proportional to y(\). HI The MFS charge points are then given by y 3 = Z(xj ~ w{Xj))- 

We sketch our motivation for the above algorithm. The curve (50) can be viewed as defining a z-plane 
curve which under some new parametrization (annular conformal map) is the image of a line Im s = C < 0. 
However, because of the curvature behavior near singularities induced by the cos terms in (50), the r for the 
new parametrization is large, of order 1. One may see intuitively that it is possible to choose a parametrization 
where any given singularity a has its imaginary part r a pushed to infinity by choosing the conformal map 
corresponding to the dQ held at electrostatic potential zero while a point charge is placed at the singularity; 

1 In practice this can easily be done numerically by solving the ODE u'(x) = Vj/(x) then spline fitting to construct an 
approximate inverse function for u. 
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Fig. 14. Interior Helmholtz BVP solution u for the analytic 5-foil domain given by the radial function r(8) = 1 + 0.3cos(5£>) 
at wavenumber k = 400, with boundary data given by v(z) = Re (z — p) _1 with p = 1 + 0.5i. Error norm is t = 2 X 10 — 11 . 
Basis size was N = 2000, and M = 2000 boundary points. Charge points are shown by dots and Schwarz branch singularities 
by '+' symbols. CPU time was 26 s to compute the coefficient vector en, and 1.05 hr to evaluate the interior solution u shown 
(8.2 X 10 J points on a grid of spacing 0.002). 

this analogy informed our choice of curve. Finally, the choice of charge point spacing approximates the result 
of choosing constant spacing on the line Ims = const under the new parametrization. This algorithm has 
been developed using intuition arising from electrostatics; we do not claim that it is the best possible choice. 
But we obtain very good results with it in our numerical experiments. The parameters we give above seem to 
work well in a wide variety of randomly-generated analytic curves (the one failure mode which the algorithm 
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does not yet guard against is self- intersection of the resulting T). 

In Fig. 12 we illustrate the performance of this method on a generalization of the crescent domain given 

by 

j rz / \ is 0.1 0.07 + 0.02z 0.2 

Generalized crescent: Zgc\ s ) — e : : 1 — : (52) 

e ls + a 5 e %s + ag e ls + a 7 

with 05 = 0.9, ct6 = —0.8 — 0.2i, 07 = —0.2 + 0.5i. The Schwarz function has three exterior pole-type 

singularities, shown by '*' symbols, with s-plane locations I/S5, 1/&6, and 1/07. The two rows of subfigures 

contrast the effect of -D max given by (51) at low vs high wavenumbers. For low wavenumber -D max is large 

and the contribution of the first term in (50) is small. This enables the curve to reach large negative Im 

s, hence large distances from <9f2 and a large point spacing, in regions away from singularities. We observe 

rapid convergence, reaching t = 10~ 14 by N — 280, and no exponential coefficient growth. As k increases, 

-Dmax drops and the first term in (50) starts to become significant (see dotted line in Fig. 12d), and has the 

effect of bringing the curve closer to the domain. As k — * 00 this term dominates and the curve tends to a 

constant distance of about 4 wavelengths from dQ everywhere on the curve. The bottom row in the figure 

shows the case k — 100. Once N = 400 (about 3.0 degrees of freedom per wavelength on the boundary), 

convergence is rapid, reaching 10 -10 at N = 525. Since the coefficient norm |a| w 10 5 convergence halts 

here. The resulting solution u is shown in Fig. 13. Values outside f2 have been included to highlight the 

manner in which the MFS points generate the field (very large coefficients are easily noticed due to the dark, 

highly oscillatory bands in these parts of T) . 

Finally we test a different shape at higher wavenumber in Fig. 14. This domain is challenging since it 

has, very close to the boundary, five exterior branch-type singularities in the Schwarz function. We also 

choose non-constant boundary data which itself has a singularity (however since it is outside V there is no 

need to include its contribution in the curve formula (50)). There are about 165 wavelengths across the 

domain diameter. An error norm of order 10 -11 is reached using an TV corresponding to 3.5 basis functions 

per wavelength on the perimeter, (since due to resonance the interior u values are of order 10 2 , but d as 

defined below (5) is of order 10 -5 ; these combine to give around 8 digits of relative accuracy in u). We note 

by contrast that for boundary element and boundary integral formulations it is commonly stated that 10 

degrees of freedom per wavelength are required for high accuracy. Away from singularities (or in domains 

with more distant singularities) we find barely more than 2 basis functions per wavelength are sufficient in 

the high k limit, similar to what was found for the disc in Section 2.4. We postpone further study of this 

limit, and the choice of -D max at high wavenumber, to future work. 



4. Conclusions 



The Method of Fundamental Solutions is a powerful tool for solving the Helmholtz BVP, but, as we have 
demonstrated, the achievable accuracy is limited by the size of the coefficient norm \a\. Therefore we have 
analysed, for the first time, the growth of \a\ with basis size N as one converges towards Helmholtz solutions 
in the disc and other analytic domains. We emphasize that it is not the growth in the condition number that 
we are studying (since this always grows exponentially with N), rather the growth in the norm of the least- 
squares coefficient vector. In the disc we have theorems on convergence rate (Thm. 3) and coefficient growth 
(Thms. 6 and 7), and for analytic domains we have corresponding conjectures (Conjs. 9 and 12) supported 
by numerical experiments in many domain shapes. These show that the success (numerical stability and 
hence high accuracy) of the MFS relies on a choice of charge curve which does not enclose any singularities 
of the analytic continuation of the solution u. These singularities are associated either with the analytic 
continuation of the boundary data, or with the Schwarz function of the domain. 

The conclusions for optimal choice of MFS charge points are as follows. For the unit disc, with concentric 
equally-spaced charge points, a radius between ^J~p and p is optimal (Remark 8), p being the radius of the 
nearest singularity in boundary data. For general analytic domains, charge points placed on a curve which 

2 Computation times are quoted for a single core of a 2 GHz Intel Core Duo laptop CPU, with 2 GB RAM, running GNU/Linux, 
coded in MATLAB. For evaluation of u the sum (2) was performed naively. The MATLAB Hankel function routine is also by 
no means optimal, requiring on average 2.3 microseconds per evaluation. 
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adapts to the singularity locations have been shown to perform very well (and vastly outperform charge 
points located using equipotential lines of the exterior conformal map). Thus knowledge of singularities 
in the Schwarz function and the boundary data are essential to achieve the best performance on analytic 
domains. 

Our experiments show that MFS is highly competitive with boundary integral equations, both in terms 
of basis size N and overall simplicity. At high wavenumber we show that in the disc asymptotically 2 basis 
functions per wavelength on the boundary are needed, and that in more complicated nonconvex domains with 
nearby singularities this need only increase to about 3.5 to achieve boundary error norms close to machine 
precision. In practice this unusually small N results in rapid solution of the basis coefficients, even at high 
wavenumber. However, as with boundary integral methods, the CPU time to evaluate the interior solution 
at roughly 10 grid points per wavelength is much larger, especially using a naive implementation of the sum 
(2) and MATLAB's Hankcl function routine. Replacing this evaluation of u with a Fast Multipole (FMM) 
summation would be a natural next step and is expected to result in a large speedup at high wavenumbers. 

We expect our findings on coefficient growth rates, and the new adaptive charge curve algorithm, to be 
easily extendable to the exterior Hclmholtz scattering problem, for which MFS has shown promise in the 
engineering community [21,11]. 
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Appendix A. Bessel function asymptotics 



We take the standard Taylor series [1] 

' z \ m X^ (-2 2 /4) fe 



k=0 y ' 

and in the large-m limit we may approximate (m + k)\ w mlm k , then recognize the power series for the 
exponential, giving 

1 / 7 \ m 2 / , 

■W*)~^(f) 1 4m , m^oo. (A.2) 

Similarly the standard series 

«.(•) = ~ (I)""" E {m -"-^ /4)k + 1 1„W2)J,„ W + OV~) (A.3) 



k=0 



with (m — k — 1)! w ml/m gives 

~ (£f "V /4m , m^oo. (A.4) 

Neither of these asymptotic forms are given in [1], however (A.2) has been recently noted in the physics 
community [22]. Combining these two in (12), using the reflection formulae, and ignoring the lower-order J 
contribution to the Hankcl function, gives (44). 

Bessel's equation u" + u'/r+ (1 — m 2 /r 2 )u — with the Liouville transformation w = r x ^ 2 u then changing 
variable to x = r/a, with a 2 = m 2 — j, gives the ODE 
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The WKBJ (or Liouvillc- Green) asymptotic approximation (Ch. 9.3 of [26]) for large parameter a is then 

(x- 2 - I)- 1 / 4 (Ae a SZ ^ x - 2 - ldx + Be- a J* ^~ ld A ; x<1 (evanescent) 

w ( x ) ~ I > f x , — r * / — \ (A.6) 

(1 - .t- 2 )- 1 / 4 ^e M Ji v^'^ + De-^Ji Vi-- 2 ^j ) ^ > x (oscillatory) 

where A, B,C,D E C are constants. Note that the integral in the evanescent region can be performed 
analytically and is — I\{x) as defined by (46); the integral in the oscillatory region is not needed since 
amplitude not phase is of interest. Since the solution w is continuous through the turning point x = 1 (even 
though (A.6) breaks down), there exist connection formulae relating the constants: 

C = e^l^A + e-™' A B, D = e l7! ' A B. (A.7) 

They can be found by comparing WKBJ to large- argument asymptotics of the Airy functions Ai and Bi 
on either side of the turning point {e.g. comparing 10.4.59 with 10.4.60, and 10.4.63 with 10.4.64, in [1], or 
using 9.3.91,92 of [26], or the more rigorous presentation of the Gans- Jeffreys formulae in Ch. 11 of [27]). 
If A = (w decaying as x decreases in the evanescent region) then \C\ = \D\ = \B\ giving an amplitude of 
2|£>|/(1 — a; -2 ) 1 / 4 in the oscillatory region. When transformed back such a solution u(r) corresponds to the 
J m (r) Bessel function. We match the asymptotic amplitude at large argument (see 9.2.1 of [1]) to 

fix B = i for all m. Hence the Bessel function has typical size 

f V-r 2 r 1/4 e Mr \ r<a 
\J m (r)\n{ 2 V (A.8) 
\(a 2 -r 2 )- l '\ r>a. 

Note that an amplitude is implied here in the oscillatory region r > a. Note also that a is defined above, 
and I a (r) < for r < a. Similarly matching the Hm\r) Hankel function large argument asymptotic gives 
|C| = 1, D = 0, so \A\ = 1 (which dominates), thus typical size 

f (a 2 - r 2 )~ V^e-Ur) r<a 
l^(r)|« ° a r 2 _ 1/4 6 (A.9) 
I (a 2 - r 2 ) 1/4 , r > a. 

These formulae have been checked against numerical evaluations of Bessel functions and accurately predict 
amplitudes or evanescent magnitudes everywhere apart from very close to the turning point r = a where 
they have a weak algebraic singularity, but still provide an upper bound. Substituting (A.8) and (A.9) into 
(12) gives the desired (45). 
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